1 Resumo

  • Período: 2018 a 2023;
  • Posto X;
  • Quantidade de observações iniciais: 365
  • Método de agrupamento: média mensal;
  • Quantidade de observações finais: 72;
  • Imputação de dados na observação dos meses listados abaixo foi feita a partir da média daquele ano;
  • Tentativa de ajuste com um modelo ARIMA;

A partir da base de dados dos preços de combustíveis na cidade de Goiânia, será feita a filtragem listada acima. Foi selecionado um único posto revendedor - o posto com maior número de observações no período citado - contendo 365 valores divididos em 188 datas distintas; a seguir, foi tomada a média de cada mês essa será a variável estudada; foi feito a imputação de valores a cinco meses que não possuíam dados (nov-2018, dez-2018, fev-2019, set-2020, fev-2023) por meio da média dos seus respectivos anos; por fim,

2 Organização dos dados

2.1 Antes

# Retirando 'Revenda' e 'Produto'.
dados = dados3[, -c(1,2)]

summary(dados)
##  Data.da.Coleta       Valor.de.Venda      mes                ano           
##  Min.   :2018-01-03   Min.   :3.699   Length:365         Length:365        
##  1st Qu.:2019-12-04   1st Qu.:4.690   Class :character   Class :character  
##  Median :2021-04-07   Median :4.990   Mode  :character   Mode  :character  
##  Mean   :2021-03-13   Mean   :5.396                                        
##  3rd Qu.:2022-09-29   3rd Qu.:5.890                                        
##  Max.   :2023-12-26   Max.   :7.870

2.2 Depois

# Meses faltantes
m <- check_missing_months(dados, "Data.da.Coleta")
## [1] 5
print(m)
## [1] "2018-11" "2018-12" "2019-02" "2020-09" "2023-02"
# xodo: 5 faltas "2018-11" "2018-12" "2019-02" "2020-09" "2023-02"
# hiper moreira: 18 faltas
# viena: 5
# Monaco: 4

Como há 5 meses sem nenhuma observação foi utilizado da imputação de dados. Para cada mês listado acima será imputado o valor da média do respectivo ano, portanto, novembro e dezembro de 2018 terão o mesmo valor.

# Calculando médias anuais
aux = dados %>% group_by(ano) %>% mutate(media_anual = mean(Valor.de.Venda))

#unique(aux[, c(4, 5)])

# "2018-11" "2018-12" "2019-02" "2020-09" "2023-02"
# Adicionar os seguintes valores
linha1 = data.frame(matrix(c("2018-11-01", 4.661141, "2018-11", "2018",
                             "2018-12-01", 4.661141, "2018-12", "2018",
                             "2019-02-01", 4.665063, "2019-02", "2019",
                             "2020-09-01", 4.460714, "2020-09", "2020",
                             "2023-02-01", 5.702877, "2023-02", "2023"),
                           byrow= T, ncol=4))
colnames(linha1) <- colnames(dados)

linha1[, 1] <- as.Date.character(linha1[, 1])
linha1[, 2] <- as.numeric(linha1[, 2])
d1 = unique(dados_imput[, c(3, 5)])

a = qcc::qcc(d1$media_mensal, type = 'xbar.one',
         title = 'Carta da média amostral (Dados finais)')

a = qcc::ewma(d1$media_mensal, title = 'Carta MMEP (Dados finais)')

a = qcc::cusum(d1$media_mensal, title = 'Carta "Soma cumulativa" (Dados finais)')

# AGRUPAMENTO TRIMESTRAL
m1 = matrix(d1$media_mensal, ncol = 3, byrow = T)
m11 = matrix(d1$mes, ncol=3, byrow = T)

a = qcc::qcc(data = m1, type = "S")

a = qcc::qcc(data = m1, type = "R")

A aplicação das cartas aos dados autocorrelacionados indica uma quantidade muito grande de pontos fora de controle. É de se esperar que o processo não seja bem medido e tais dados façam com que as cartas tornem-se, frequentemente, errôneas.

3 Série temporal

par(bg='gray')
d1 = unique(dados_imput[, c(3, 5)])

serie = ts(d1[, 2], start=c(2018, 1), frequency = 12)
#serie = ts(dados_imput$media_mensal, start=c(2018, 1), end=c(2023, 12), frequency = 12)
plot(serie, ylab="Valor", xlab="Obs.")

forecast::autoplot(decompose(serie))

Verificando a presença de autocorrelação entre as observações da série nota-se que há o decaimento exponencial à medida que se aumenta o lag.

par(bg='gray')
acf(d1$media_mensal, main = "Correlação entre as médias mensais")

4 Modelo ARIMA

Após toda organização do conjunto de dados, foram testados diversos modelos SARIMA. A seguir há uma relação dos modelos que apresentaram resíduos com: normalidade, estacionariedade e independência.

Outros modelos que também apresentaram todos pressupostos.

#aux = forecast::auto.arima(serie)

ar_ = forecast::Arima(serie, order = c(2, 3, 0),
                      seasonal = c(0,0,0), method = "ML")

summary(ar_)
## Series: serie 
## ARIMA(2,3,0) 
## 
## Coefficients:
##           ar1      ar2
##       -0.9819  -0.4259
## s.e.   0.1076   0.1073
## 
## sigma^2 = 0.2455:  log likelihood = -48.96
## AIC=103.93   AICc=104.3   BIC=110.63
## 
## Training set error measures:
##                       ME      RMSE      MAE       MPE    MAPE      MASE
## Training set 0.001505884 0.4779935 0.353291 0.1276067 6.67558 0.3425822
##                     ACF1
## Training set -0.09204694
# AMBOS PARÂMETROS DO MODELO ARIMA SÃO NEGATIVOS!!

shapiro.test(residuals(ar_))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(ar_)
## W = 0.96813, p-value = 0.06461
tseries::adf.test(residuals(ar_))
## Warning in tseries::adf.test(residuals(ar_)): p-value smaller than printed
## p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  residuals(ar_)
## Dickey-Fuller = -6.1038, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
Box.test(residuals(ar_), type="Ljung")
## 
##  Box-Ljung test
## 
## data:  residuals(ar_)
## X-squared = 0.63581, df = 1, p-value = 0.4252
forecast::checkresiduals(ar_)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,3,0)
## Q* = 23.918, df = 12, p-value = 0.02087
## 
## Model df: 2.   Total lags used: 14

5 Aplicação das cartas de controle

#agrup = qcc::qcc.groups(as.vector(ar_$residuals),
#                        as.vector(unique(dados_imput[, 3])))

g1 = qcc::qcc(ar_$residuals, type = "xbar.one", plot = T,
              xlab="Grupo (Mês)",
              ylab="Estatística do grupo",
              add.stats = F,
              label.limits = c("LCI", "LCS"))

g1$violations$beyond.limits
## [1] 57
g2 = qcc::ewma(ar_$residuals)

g3 = qcc::cusum(ar_$residuals)

#plot(g1, restore.par = F)
#head(d1)
d1[order(d1$media_mensal, decreasing = T), ]
## # A tibble: 72 × 2
## # Groups:   mes [72]
##    mes     media_mensal
##    <chr>          <dbl>
##  1 2022-05         7.63
##  2 2022-04         7.62
##  3 2022-06         7.61
##  4 2022-03         7.53
##  5 2022-01         7.36
##  6 2021-11         7.34
##  7 2021-12         7.03
##  8 2021-10         7.02
##  9 2022-02         6.99
## 10 2021-09         6.66
## # ℹ 62 more rows
d1[56:58, ]
## # A tibble: 3 × 2
## # Groups:   mes [3]
##   mes     media_mensal
##   <chr>          <dbl>
## 1 2022-08         5.10
## 2 2022-09         4.96
## 3 2022-10         4.99
g1$statistics[56:58]
##        56        57        58 
## 0.9355005 1.6073610 0.3819759
g1 = qcc::qcc(ar_$residuals, type = "xbar.one", plot = F,
              add.stats = F)

# 1. Convert qcc object to a data frame for ggplot2
qcc_data <- data.frame(
  x = 1:length(g1$data), # x-axis (sample number)
  Value = g1$data,      # y-axis (data values)
  CL = mean(g1$data[, 1]),     # Center Line
  UCL = g1$limits[2],    # Upper Control Limit
  LCL = g1$limits[1],    # Lower Control Limit
  violations = ifelse(g1$data > g1$limits[2] | g1$data < (g1$limits[1] - 3*(g1$limits[2] - g1$limits[1])), "Violação", "Sob controle"))

# 2. Create the ggplot2 plot
gg1 <- ggplot(qcc_data, aes(x = x, y = Value)) +
  theme_classic() +
  geom_point(aes(color = violations)) +  # Points, colored by violations
  geom_line() +                      # Connect the points
  geom_hline(yintercept = qcc_data$CL[1], linetype = "solid", color = "black",
             linewidth = .1) + # Center line
  geom_hline(yintercept = qcc_data$UCL[1], linetype = "dashed", color = "red",
             linewidth = .3) + # Upper control limit
  geom_hline(yintercept = qcc_data$LCL[1], linetype = "dashed", color = "blue",
             linewidth = .3) + # Lower control limit
  labs(title = "", x = "Grupo (Mês)", y = "Estatística do grupo") + # Labels
  scale_color_manual(values = c("Violação" = "red", "Sob controle" = "black")) + # Custom colors
  annotate(geom = 'text', label = c(paste0("LSC = ", round(qcc_data$UCL[1], 3)),
                                    paste0("LC = ", round(qcc_data$CL[1], 3)),
                                    paste0("LIC = ", round(qcc_data$LCL[1], 3))),
           size = 3, x = c(5,1,5), y = c(qcc_data$UCL[1] + .1 ,
                                         qcc_data$CL[1] + .1,
                                         qcc_data$LCL[1] + .1)) +
  theme(legend.position = "none") # Remove the legend

plotly::ggplotly(gg1)

5.1 Adicionais

Gráficos utilizados no relatório.

ar1 = matrix(ar_$residuals, ncol = 3, byrow = T)

ar11 = matrix(d1$mes, ncol = 3, byrow = T)
ar111 = matrix(c(1:72), ncol = 3, byrow = T)

am1 = qcc::qcc(ar1, type = "R")

g1$violations$beyond.limits
## [1] 57
#ar11[g1$violations$beyond.limits, ]

#qcc::qcc(ar1, type="S")
ar1 = matrix(ar_$residuals, ncol = 4, byrow = T)

ar11 = matrix(d1$mes, ncol = 4, byrow = T)
ar111 = matrix(c(1:72), ncol = 4, byrow = T)

am1 = qcc::qcc(ar1, type = "R")

am1$violations$beyond.limits
## [1] 14 15
#ar11[am1$violations$beyond.limits, ]

qcc::qcc(ar1, type="S")

## List of 11
##  $ call      : language qcc::qcc(data = ar1, type = "S")
##  $ type      : chr "S"
##  $ data.name : chr "ar1"
##  $ data      : num [1:18, 1:4] 0.00104 0.10152 0.28529 0.16085 0.04526 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics: Named num [1:18] 0.0244 0.4928 0.363 0.2921 0.2695 ...
##   ..- attr(*, "names")= chr [1:18] "1" "2" "3" "4" ...
##  $ sizes     : int [1:18] 4 4 4 4 4 4 4 4 4 4 ...
##  $ center    : num 0.457
##  $ std.dev   : num 0.496
##  $ nsigmas   : num 3
##  $ limits    : num [1, 1:2] 0 1.04
##   ..- attr(*, "dimnames")=List of 2
##  $ violations:List of 2
##  - attr(*, "class")= chr "qcc"
ar1 = matrix(ar_$residuals, ncol = 3, byrow = T)

ar11 = matrix(d1$mes, ncol = 3, byrow = T)
ar111 = matrix(c(1:72), ncol = 3, byrow = T)

am1 = qcc::qcc(ar1, type = "S")

g1$violations$beyond.limits
## [1] 57
#ar11[g1$violations$beyond.limits, ]
# 1. Convert qcc object to a data frame for ggplot2
qcc_data <- data.frame(
  x = 1:length(am1$statistics), # x-axis (sample number)
  Value = am1$statistics,      # y-axis (data values)
  CL = mean(am1$statistics),     # Center Line
  UCL = am1$limits[2],    # Upper Control Limit
  LCL = am1$limits[1],    # Lower Control Limit
  violations = ifelse(am1$statistics > am1$limits[2] | am1$statistics < (am1$limits[1] - 3*(am1$limits[2] - am1$limits[1])), "Violação", "Sob controle"))

# 2. Create the ggplot2 plot
gg1 <- ggplot(qcc_data, aes(x = x, y = Value)) +
  theme_classic() +
  geom_point(aes(color = violations)) +  # Points, colored by violations
  geom_line() +                      # Connect the points
  geom_hline(yintercept = qcc_data$CL[1], linetype = "solid", color = "black",
             linewidth = .1) + # Center line
  geom_hline(yintercept = qcc_data$UCL[1], linetype = "dashed", color = "red",
             linewidth = .3) + # Upper control limit
  geom_hline(yintercept = qcc_data$LCL[1], linetype = "dashed", color = "blue",
             linewidth = .3) + # Lower control limit
  labs(title = "", x = "Grupo (Mês)", y = "Estatística do grupo") + # Labels
  scale_color_manual(values = c("Violação" = "red", "Sob controle" = "black")) + # Custom colors
  annotate(geom = 'text', label = c(paste0("LSC = ", round(qcc_data$UCL[1], 3)),
                                    paste0("LC = ", round(qcc_data$CL[1], 3)),
                                    paste0("LIC = ", round(qcc_data$LCL[1], 3))),
           size = 3, x = c(5,1,5), y = c(qcc_data$UCL[1] + .1 ,
                                         qcc_data$CL[1] + .1,
                                         qcc_data$LCL[1] + .1)) +
  theme(legend.position = "none") # Remove the legend

plotly::ggplotly(gg1)